Direct multiscale coupling of a transport code to gyrokinetic turbulence codes 
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Direct coupling between a transport solver and local, nonlinear gyrokinetic calculations using the 
multiscale gyrokinetic code TRINITY [M. Barnes, Ph.D. thesis, arXiv:0901.2868 is described. The 
coupling of the microscopic and macroscopic physics is done within the framework of multiscale 
gyrokinetic theory, of which we present the assumptions and key results. An assumption of scale 
separation in space and time allows for the simulation of turbulence in small regions of the space-time 
grid, which are embedded in a coarse grid on which the transport equations are implicitly evolved. 
This leads to a reduction in computational expense of several orders of magnitude, making first- 
principles simulations of the full fusion device volume over the confinement time feasible on current 
computing resources. Numerical results from TRINITY simulations are presented and compared with 
experimental data from JET and ASDEX Upgrade plasmas. 



I. INTRODUCTION 

A fundamental challenge of fusion science is to max- 
imize fusion power, which is determined primarily by 
macroscopic profiles of density and temperature. These 
profiles, which vary spatially on the system scale and 
evolve on the energy confinement time scale, drive turbu- 
lence at micro-scales in space and time. In the absence 
of MHD instability, this micro-turbulence is the domi- 
nant source of heat flux observed in standard tokamaks, 
which sets rigid constraints on the macroscopic profiles. 1 
Consequently, it is of critical importance to understand 
the self-consistent interaction between the macroscopic 
profiles and the micro-turbulence. 

This is a challenging task due to both the wide range of 
scales involved and the high dimensionality of the system. 
The electron turbulence space scale, which is comparable 
to the electron Larmor radius, will be on the order of 0.1 
millimeters for a device such as ITER, 2 -' which has a mi- 
nor radius of about two meters. Similarly, the electron 
turbulence time scale is approximately one microsecond, 
much smaller than the expected energy confinement time 
of several seconds. Additionally, the instabilities driv- 
ing the turbulence are kinetic in nature, requiring treat- 
ment of the velocity space in addition to the configuration 
space. Including all of these dynamics directly in a single 
simulation is not feasible on current computing resources. 

As long as all relevant dynamics occur on time scales 
long compared to particle gyro-motion, it is possible to 
average over the gyro-orbits and eliminate the gyro-angle 
as a phase space variable. Furthermore, for magnetic con- 
finement devices with sufficiently small (ratio of ion 
Larmor radius to plasma minor radius) , there is expected 
to be a separation between micro- and m acro-s cales in 
space and timeP^The ^/-gyrokinetic modeP^ exploits 
these scale separations, simplifying the system consider- 
ably. Given a set of fixed macroscopic profiles, it allows 
for the calculation of turbulent fluxes in a 5-D phase 
space. Such calculations have been performed numeri- 



cally in gyrokinetic codes for more than two decades, pro- 
viding much insight into the nature of kinetic instabilities 
and micro-turbulence. First-principles ^/-gyrokinetic 
simulations have steadily advanced in their sophistication 
and physical fidelity and have become routine in recent 
years. 

However, these codes only model the effect of macro- 
scopic profiles on micro- turbulence: They do not provide 
quantitative information on how this turbulence subse- 
quently affects the evolution of the macroscopic profiles. 
Accurately calculating fluxes from experimental profiles 
is also problematic because of the acute sensitivity of 
the fluxes to small changes in the input profile gradients 
(which would arise due to uncertainty in experimental 
measurements) P In an attempt to address these issues, 
full-/ gyrokinetic codes have been developed, which do 
not explicitly assume the scale separations listed above. 
Subsequently, they are able to take into account the two- 
way interaction between turbulence and equilibrium ther- 
modynamic profiles. However, well-resolved full-/ sim- 
ulations would be extremely expensive numerically be- 
cause of the wide range of scales described above and 
because of the necessity of calculating the distribution 
function and fields to very high order. It has also re- 
cently been argued that this approach as currently for- 
mulated leads to an unphysical source of toroidal angular 
momentumP 

An alternative approach commonly used to study the 
interaction of the micro- and macro-physics is to solve 
fluid transport equations with a reduced model for the 
turbulent fluxes. Typically, the transport is modeled 
as a diffusive process, with turbulent diffusivities com- 
ing from a wide range of models, including empirical 
fits to experiment o r simulation and theory-based esti- 
mates.^ 10 * 11 * 12 * 13 ^ These reduced models have provided 
a basic qualitative understanding of the multiscale in- 
teraction and are capable in some cases of giving good 
quantitative agreement with first-principles, nonlinear 
gyrokinetic calculations! 15 * 16 ' However, such models do 
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not permit detailed validation studies because they do 
not produce fluctuation spectra or other data that can 
be experimentally challenged. Without careful valida- 
tion, even first-principles reduced models (with no ad- 
justable parameters) cannot be fully trusted for predic- 
tions of performance in new operational regimes, and re- 
duced models with fit parameters adjusted for current 
experiments have even less credibility. There are known 
cases for which even first-princ iples reduced models cur- 
rently available are inadequate! 17 1 18 1 

Consequently, it is desirable to couple not only re- 
duced flux models, but also direct numerical simulations 
of the turbulence, to transport solvers. This has been 
done in Ref. [19] for a Hasegawa-Wakatani 2D fluid model 
of the turbulence using an implicit relaxation technique 
and allowing for nonlocality. Here, we describe coupling 
to local, nonlinear, 5D gyrokinetic turbulence calcula- 
tions using a Newton method (similar to that described 
in Ref. |2Q|) , which accelerates convergence by more than 
an order of magnitude for typical parameters. This cou- 
pling is achieved using the multiscale gyrokinetic code, 
TRINITYj 2 -^ which can use nonlinear fluxes from the con- 
tinuum gyrokinetic codes GS^P^ or GENE. 23 24 Our ap- 
proach is similar to that employed in TGYROp^ with the 
key distinction that TRINITY evolves the macroscopic 
profiles in time, whereas TGYRO assumes a steady-state 
and solves the volume-integrated transport equations for 
profile gradients. 

It should be noted that some meso-scale phenomena in 
space and time are not formally considered in the stan- 
dard multiscale gyrokinetic model. Our simulations ig- 
nore low-order magnetic islands and so are directly ap- 
plicable only during MHD-quiescent periods of plasmas. 
Rapid cold/heat pulse propagation (such as following a 
sawtooth or ELM crash) is possible in the TRINITY code 
because of the presence of stiff critical gradients in ITG 
and TEM turbulence, though the transport time step 
would of course have to be reduced in TRINITY during 
such a transient event to be able to follow its propa- 
gation. It should be noted that flux-tube simulations 
include the contribution to the heat flux of avalanches 
on all scales up to the radial size of the flux tubes. If 
even longer wavelength avalanches were important, then 
the heat flux would increase as the flux-tube simulation 
domain was made larger, so convergence studies can be 
used to test this. PIC and continuum flux-tube simu- 
lations generally find that the flux converges with suf- 
ficiently large simulation size, of order the sizes we are 
using here. Previous gyrokinetic studies have found that 
some modest non-local turbulence spreading may occur 
over distances of a few radial eddy sizes™ but in the 
core region of the large tokamaks we are studying here 
(and even more so at reactor scales) this should usually 
be a small effect. It should be acknowledged that the 
separation of scales assumed for the core plasma in this 
paper may break down in the edge region of the plasma 
because gradient scale lengths and eddy sizes may not be 
very different near the edge, so non-local effects may be 



important there. 

This paper is organized as follows: in the next sec- 
tion we state the fundamental assumptions of the multi- 
scale model and present the closed system of equations 
that results from gyrokinetic expansion of the Maxwell- 
Boltzmann system. In Sec. |III[ we describe the numer- 
ical scheme used in TRINITY to simulate the multiscale 
gyrokinetic system of equations. We also give estimates 
for the space and time domain savings provided by the 
multiscale scheme. Sec. [IV] contains results from TRINITY 
simulations of L-mode and H-mode discharges from JET 
and ASDEX Upgrade. We show that the numerical data 
from TRINITY is in good quantitative agreement with re- 
constructions of experimental data. Finally, we conclude 
in Sec. [V] with a summary and a discussion of possible 
future directions for research. 



II. THEORETICAL FRAMEWORK 

In this section, we state the fundamental assumptions 
present in our multiscale model and present the resulting 
closed system of equations that must be solved. These 
equations, which are a rigorous asymptotic limit of the 
full Maxwell-Boltzmann system, have been derived in de- 
tail in Refs. [271 and [28| We include a brief overview of 
the key results here for completeness. 

As the starting point for our analysis, we begin with 
the coupled system consisting of Maxwell's equations and 
the driven Fokker-Planck equation: 

^ = C[f s ] + S s [f s ], (1) 

where f s = / s (x, v,£) represents the distribution of par- 
ticles of species s in position (x) and velocity (v) space, 
C represents the effect of two-particle Coulomb interac- 
tions, and S s is a source term accounting for the external 
injection of particles, momentum, and energy. This sys- 
tem of equations describes all of the important dynamics 
in fusion plasmas and is consequently intractable, both 
analytically and numerically. Since we are interested in 
studying the interaction of the plasma micro-turbulence 
with the macroscopic profiles, we simplify the system 
by adopting a variant of the standard ^/-gyrokinetic or- 
dering. 6 This ordering imposes constraints on the rela- 
tive amplitudes and space-time scales of the macro- and 
micro-physics. 

We first decompose the distribution function into 
macroscopic and microscopic quantities by ensemble av- 
eraging: / = (/) +5/, with the angled brackets denoting 
an ensemble average. Defining the smallness parameter 
to be e = p/L, where p is the Larmor radius and L is 
a macroscopic scale length, we order each of the terms 
in the Maxwell-Boltzmann equations. The assumptions 
employed in ordering the terms are as follows: (1) The 
fluctuations are assumed to be low amplitude compared 
to macroscopic quantities, such that Sf ~ e(f). This 
is in good agreement with core measurements from a 
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number of modern fusion experiments, which find den- 
sity and temperature fluctuation s < 1% of the macro- 
scopic densities and temperatures (2) The micro- 
turbulence is assumed to be spatially anisotropic with 
macro-scale variations along and micro-scale (i.e. Larmor 
radius) variations across the equilibrium magnetic field. 
Experimental measurements of turbulence parallel and 
perpendicular correlation lengths support this hypothe- 
s - g j 29 | 32 | 33 | ^ All frequencies of interest are assumed to 
be well below the ion cyclotron frequency, and the evolu- 
tion of the macroscopic profiles is taken to be much slower 
than the turbulent fluctuations (d(f)/dt ~ edSf/dt). 
Again, this transport ordering is in agreement with ex- 
perimental evidence^ (4) We order d5f/dt ~ C[Sf], 
with the collision frequency, v, defined such that v > 
ed\nSf/dt. Consequently, Sf is allowed characteristic 
scales in velocity space, Sv, of size y/evth ^ Sv < v t h- 
This collision frequency ordering is satisfied even in very 
collisionless plasmas such as anticipated in ITER; (5) 
Macroscopic flows are assumed to be comparable to the 
ion thermal speed, with microscopic flows (i.e. E x B 
velocity) taken to be much smaller. For simplicity, we 
consider in this paper the case in which the Mach num- 
ber of the macroscopic flow is taken to be small as a 
subsidiary ordering; (6) The external particle, momen- 
tum, and energy sources are assumed to affect the system 
evolution on the confinement time scale, consistent with 
experiment. 

Expanding the distribution function and fields in e and 
applying the above ordering assumptions to the Maxwell- 
Boltzmann system results in a hierarchy of equations 
that is ultimately closed by ensemble and flux surface 
averaging and taking moments of the evolution equation 
for the lowest order (macroscopic) distribution function, 
fo — (fo)- One finds that A is a gyroangle-independent, 
shifted Maxwellian whose evolution is governed by the 
following transport equations: 
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= -(H s ) + -n s J2^u(Tu-T s ) + (S p ), (4) 

where R is the major radius, <p is the physical toroidal 
angle, ip is the flux label, the overline denotes a flux sur- 
face average, and V = dV/dip, with V being the volume 
enclosed by the flux surface. The evolved quantities n s , 
p s , and L are the ensemble- averaged density, pressure, 
and species-summed toroidal angular momentum, respec- 
tively. In the low Mach limit we are considering, the den- 
sity and pressure are constant on flux surfaces. Terms 



denoted by S represent external sources, with subscripts 
indicating the relevant injected quantity. The collisional 
energy exchange frequency, v £ su is given in Ref. 34, The 
terms T, 7r, Q, and H are fluxes and heating generally 
consisting of classical, neoclassical, and turbulent contri- 
butions. They are given by 

T = V^-Jd 3 v (v x 5h + v B (A) + pC [p • V/o]) (5) 
7T = W • J d 3 v (mR 2 v • V0) v x Sf t (6) 
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where x — ^ ~ v ' ^A/c is the generalized electro- 
magnetic potential fluctuation, v x = c/Bb x Vx is 
the particle drift due to the fluctuating fields, v# = 

(b/fi) x + vjjb • Vb) contains the magnetic 

drifts, e is the species charge, p is the gyroradius vec- 
tor, and D/Dt = d/dt + u • V, with u = Rw(^)(j) the 
equilibrium flow velocity. 

The components of the first order distribution func- 
tion, A, and the fluctuating potentials are obtained by 
solving the neoclassical and gyrokinetic equations, cou- 
pled to the low-frequency Maxwell's equations. These 
equations are all self-consistently obtained as part of the 
multiscale gyrokinetic expansion. The neoclassical equa- 
tion governing (A) is 

——v\\b • — — . (9 
cT 11 dt v J 

The gyrokinetic equation determining the evolution of 

Sfi is 

^- + (v\\h + (v x ) R + w B ) ■ Vft - (C[h]) n 

fofD(x) n ^/ v Id" din A 



^[(A)]-^i|b-v(A)=v B -VA- 
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where 8h = h- (q6$/T)f , I{$) = q(ip)/(l/R 2 ) is the 
toroidal flux function, q(ip) is the safety factor, (.) R de- 
notes a gyroaverage at fixed guiding center position, R, 
and the subscript on the D/Dt operator indicates that 
it is to be evaluated at R. The low-frequency Maxwell's 
equations are given by 
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The final elements needed to close this system are 
equations for the evolution of the magnetic geometry. In 
particular, the magnetic field can be specified in an ax- 
isymmetric system if given the poloidal flux, ip, and the 
toroidal flux function, I(ip). The toroidal flux function 
is evolved by taking the toroidal component of Faraday's 
Law and flux surface averaging: 
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with the B • 9Aq / dt term obtainable from the neoclas- 
sical equation (Eq. [9|. Eqs. [9] and [l4] are then coupled 
to the Grad-Shafranov equation, which uses the updated 
macroscopic pressure from Eq. [4] to obtain ip and close 
the system: 



R 2 V 



(Y±\ = 



l R 2 



dip 



4 ^ 2 E 



dp s 
dtp 



(15) 



III. NUMERICAL METHOD 

We describe in this section the numerical model we 
have developed for solving the system described by 
Eqs. [2]|T5j This model is implemented in the multiscale 
gyrokinetic transport solver TRINITY.^ Currently, there 
are a few additional assumptions employed in TRINITY 
to simplify the multiscale system presented in the pre- 
vious section. In what follows, we provide a numerical 
prescription for solving the full system, pointing out the 
places where the TRINITY model has been simplified. 

A simple sketch of our multiscale numerical model is 
given in Fig. [I] A direct numerical simulation would 
require us to have a fine space-time mesh over the full 
device volume and over at least a confinement time. We 
exploit scale separation present in the system to dras- 
tically reduce the domain over which a fine mesh is re- 
quired. Our assumption of time scale separation between 
the turbulence and the equilibrium allows us to fix equi- 
librium quantities while we evolve the turbulence to sat- 
uration. Additionally, it allows us to use steady-state, 
time-averaged fluxes in our transport equations. Conse- 
quently, we need only resolve turbulence time scales for 
short periods of time, between which we can take large 
time steps characteristic of the confinement time. As the 
separation of scales gets wider, the simulation domain 
savings from this approach grows: The cost of simulat- 
ing small p* devices is no greater than that for moderate 
devices. The time domain savings for a device like 
ITER is a factor of hundreds. 

Similarly, spatial scale separation allows us to assume 
that macroscopic quantities (and their associated gradi- 
ent scale lengths) are constant across the radial domain 
in which we simulate turbulent dynamics. As long as 
the turbulence simulation domain is wide enough in each 
dimension, the turbulence at the ends of the domain is 
uncorrelated. Statistically periodic boundary conditions 



then apply. The result of this local approximation is a 
flux tube simulation domain for the turbulence (Fig. [I]), 
which can be used to periodically map out a flux sur- 
face. Comparisons between local and global gyrokinetic 
simulations have shown that the local approximation is 
valid for small p*, 35 as it must be for the gyro-Bohm 
scaling suggested by high confinement experiments 3 -^ to 
hold. Once again, the spatial domain savings increases 
with the scale separation. On small devices, the sim- 
ulation volume is reduced by approximately a factor of 
one hundred. 



A. Discretization of the transport equations 



The transport equations (Eqs. |2]|4[ ) are stiff, nonlinear 
partial differential equations. In order to take the large 
time steps required by our multiscale scheme, we must 
treat them implicitly. We allow for a general, single-step 
time discretization, but we primarily use first-order back- 
wards differences for steady-state systems and second- 
order backwards differences for time-dependent systems. 
An adaptive time step is employed, allowing for accu- 
rate time evolution with large time steps. The nonlinear 
terms are treated implicitly by linearizing them using 
a standard, multi-iteration Newton's method similar to 
that given in Ref. [20j For instance, the normalized heat 
fluxes at the m + 1 time level, which are nonlinear func- 
tions of the macroscopic profiles, are expanded as follows: 
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where p denotes the iteration index within each time step, 
Q = (Q/pv t h)(a/ p) 2 j and all quantities are understood 
to be evaluated at time level m + 1. The vector y con- 
tains the profiles of the fundamental macroscopic, time- 
dependent quantities in the simulation. This consists 
of the two free flux functions from the Grad-Shafranov 
equation, ip and I(ip), as well as the species density and 
pressure and the species-summed toroidal angular mo- 
mentum. We are not currently evolving the magnetic 
equilibrium in TRINITY, so that ip and I(ip) are fixed in 
time. 



Discretizing Eq. 16 in space, we obtain 
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where the subscript denotes the spatial index. In the 
local approximation, the fluxes depend only on the local 
values of macroscopic quantities and their gradients. The 
above expression thus reduces to 
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FIG. 1: (Left:) Cartoon of multiscale space-time grid used in TRINITY. Blue (vertical) and yellow (horizontal) regions represent 
the radial and time domains, respectively, over which a fine mesh is used to calculate turbulent fluxes. The overlapping regions, 
colored in green, denote the reduced space-time domain used in TRINITY. These green patches, each of which corresponds to 
a nonlinear gyrokinetic flux tube simulation (right), are grid points in the coarse space-time mesh used to solve the transport 
equations. 



In TRINITY, we make the further simplifying assumption 
that the fluxes depend more strongly on profile gradi- 
ents than the local values themselves. We then neglect 
dQ/dy, giving us the final expression 
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(19) 



If this assumption is not satisfied, it affects only the rate 
of convergence of the solution, not its accuracy. 

In order to numerically calculate dQj/dy'^ we must 
employ finite differences. This requires us to compute 
Qj at multiple values of yj. This is equivalent to cal- 
culating the fluxes for both the nominal profiles and for 
additional profiles corresponding to each gradient that 
must be perturbed. The total number of flux tube cal- 
culations required during each transport time step, TV, is 
given by 



TV = n r x (1 + n p ), 



(20) 



where n r is the number of radial grid points in the trans- 
port solver and n p is the number of macroscopic profiles 
being evolved. 

Once the transport equations are linearized, it is 
straightforward to implicitly evolve them. The expense 
of the implicit evolution is negligible when compared 
with the cost of the nonlinear turbulence calculations so 
that there is no need to use an approximate Jacobian. 
Because inversion of the Jacobian is essentially free, a 
dense matrix, arising from high order spatial derivatives, 
is tractable. 



B. Schematic 

To begin the numerical calculation, the initial state 
of the plasma must be specified. In particular, enough 
information must be given to calculate the fluxes and 
heating from Eqs. [5](8] This requires local information 
about the magnetic equilibrium, as well as values for the 
macroscopic density, flow, and temperature and their gra- 
dients at each of the flux surfaces comprising the radial 
grid for the transport solver. TRINITY is currently ca- 
pable of both an analytic and numerical specification of 
these quantities, with experimental values taken from the 
publicly accessible ITER profile database.^ Once these 
quantities are obtained at each of the radial grid points, 
TRINITY calls a solver for the fluxes. For the ion neo- 
classical heat flux, TRINITY currently uses the simplified 
analytic model given in Ref. [37j All other neoclassical 
and classical fluxes, which are typically small compared 
to turbulent fluxes J3HCSD are neglected. For a more accu- 
rate treatment of neoclassical effects, one could interface 
to a code such as that given in Ref. HOI which solves Eq.|9] 
directly. This will be the subject of future work. 

For the turbulent fluxes, there are interfaces within 
TRINITY to two widely - used, nonlinear gyrokinetic codes, 
GS^ and GENE pMl Additionally, there are options to 
use the IFS-PPPL modeP' and other simpler analytic 
models. Because we are using the local approximation, 
the flux calculations at each radius are independent for a 
given transport time step. Consequently, each flux tube 
calculation can be run in parallel, with the only com- 
munication occurring when the fluxes are gathered to 
advance the transport equations. Since each flux tube 
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TABLE I: Experimental parameters 



shot 


time (s) 


B* (T) 


a (m) 


I P (MA) 


JET 19649 


8.7 


3.12 


1.16 


3.05 


JET 42982 


14.8 


4.0 


0.95 


3.76 


AUG 19649 


1.35 


2.48 


0.48 


1.0 



calculation parallelizes with high efficiency to thousands 
of processors, our scheme with tens of flux tubes can eas- 
ily scale to hundreds of thousands of processors. 

Once the steady-state turbulent fluxes are computed, 
they are time- and flux tube- averaged and passed back 
to the transport solver. The discretized transport equa- 
tions are then solved to obtain updated densities, pres- 
sures, and the species-summed toroidal angular momen- 
tum. Boundary conditions are required to obtain unique 
solutions. At the outermost radius in the simulation, we 
fix the values of the thermodynamic profiles, typically 
to values taken from experiment. As the other bound- 
ary condition, we take the product of V with the fluxes 
to vanish at the magnetic axis. Within each transport 
time step, the equations are iterated until the relative er- 
ror upon successive iterations is less than a user-specified 
tolerance. In practice, we find that two iterations is suf- 
ficient to obtain accurate results. 

Currently, we only consider static magnetic equilibria 
in TRINITY. This eliminates several terms in Eqs.[2]|4j and 
foregoes the necessity of solving Eqs.[l4]and[l5] This ap- 
proximation is strictly valid in the limit of (3 = Snp/B 2 ^> 
^m e /mi, where the magnetic geometry evolves on a re- 
sistive time which is much longer than the energy con- 
finement time. However, if one were to evolve the mag- 
netic equilibrium, then the next step would be to use the 
parallel current obtained from Eq. [9] in Faraday's Law 
(Eq. 14) to calculate When combined with the up- 

dated pressure gradient, this allows us to solve for ip in 
With the updated thermodynamic profiles and 



Eq. 15 



magnetic equilibrium, we complete the feedback loop by 
solving for updated fluxes. This process is repeated as 
many times as is necessary to evolve the macroscopic pro- 
files beyond the time of interest. For steady-state simu- 
lations, approximately ten to fifteen transport time steps 
are typically required in TRINITY for convergence. 



IV. SIMULATION RESULTS 

We illustrate the utility of our multiscale model in 
this section by presenting numerical results from TRINITY 
simulations and comparing these results with experimen- 
tal data. We consider an L-mode discharge from JET 
(shot #19649) and H-mode discharges from JET (shot 
#42982) and ASDEX Upgrade (shot #13151). For sim- 
plicity, we consider time slices taken from approximately 
steady-state periods of each discharge. Initial thermo- 
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FIG. 2: Comparison of steady-state density and tempera- 
ture profiles constructed from JET shot #19649 by TRANSP 
(points and dotted lines) with those calculated in TRINITY 
(solid lines). 



dynamic profiles and external sources used in TRINITY 
are taken during these steady-state time slices from the 
TRANSfED or ASTRiP^l reconstructions provided in the 
ITER profile database. Some key experimental parame- 
ters for these shots are given in Table [IJ In all gyrokinetic 
simulations, we consider electrostatic turbulence with gy- 
rokinetic ions and electrons. 

First, we consider JET shot #19649P This was a 
standard JET L-mode pulse with 9.2 MW of neutral 
beam heating. For the TRINITY simulation, nonlinear 
fluxes from GS2 were used in the transport solver. We 
used a Miller local equilibrium model 44 for the magnetic 
geometry in these GS2 simulations, with the necessary 
parameters taken from the ITER profile database. In 
physical space, we used 16 grid points along the equilib- 
rium magnetic field and a 40 x 25 grid in the perpen- 
dicular plane, with perpendicular box widths at the out- 
board midplane of 64 pi. In the dealiased Fourier space, 
this corresponds to covering \k$Pi\ =0,0.1,0.2, ...0.8 and 
\k r pi\ = 0, 0.1, 0.2, 1.3. The velocity space grid con- 
sisted of 10 velocities and 32 pitch angles, giving ap- 
proximately 10 7 mesh points for each two-species GS2 
simulation. These simulations employed a hyperdiffusion 
operator whose magnitude is scaled by the shearing rate 
of the turbulence, which provides a sub-grid model of 
the cascade of fluctuations to smaller scales!™™ Previ- 
ous tests of these kinds of sub-grid models have found 
that they can reduce the needed resolution, but more de- 
tailed convergence studies in the future could be useful. 
Note also that because of the nonlinearities and stiffness 
of the transport, changes in the turbulence level of a few 
tens of percent at fixed temperature gradient have little 
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FIG. 3: Comparison of experimental (dotted lines) and simu- 
lated (solid lines) gradient scale lengths for JET shot #19649. 



effect on the final self-consistent temperature profile. 

Each flux tube simulation was run for 10 4 time steps, 
corresponding to average physical simulation times at 
each radius of approximately 400-1400 L^A^i- Elec- 
tron density and ion and electron pressures were evolved 
at 8 radial locations, giving a total of 32 flux tube sim- 
ulations per transport time step. The total number of 
mesh points required for each transport time step was 
thus approximately 3 x 10 8 . With a total of 15 transport 
time steps taken, the simulation lasted approximately 4 
hours on 5760 CRAY XT4 processors. 

A comparison of the steady-state profiles calculated 
in TRINITY with the TRANSP-reconstructed experimen- 
tal profiles is given in Fig. [2] We find good agree- 
ment for all profiles across the simulation domain (r/a = 
0.053 — 0.8), with an RMS relative error amongst all pro- 
files (^/(^ + +<4J/3) of 12% (Table Eli. The total 

and incremental stored energies differ from TRANSP val- 
ues by 9% and 12%, respectively. A comparison of the 
gradient scale lengths is given in Fig. [3} For r/a < 0.6, 
the ion temperature gradient scale lengths from TRINITY 
and TRANSP match almost perfectly. However, the elec- 
tron density and temperature gradient scale lengths do 
not agree as well, despite reasonable profile agreement. 
This illustrates a difficulty in using experimental pro- 
file measurements in standalone gyrokinetic turbulence 
calculations: small changes in experimental profiles can 
lead to significant changes in experimental gradient scale 
lengths, which drastically affect the turbulent fluxes. An 
example of this is seen in Fig. |4j Here, we compare the 
volume integrated source terms to the flux surface inte- 
grated fluxes, which should be equal in steady state. We 
see that this balance is satisfied for the self-consistent 



FIG. 4: Power balance for JET shot #19649. Solid lines are 
the steady- state, flux surface integrated fluxes calculated in 
GS2 at the end of the TRINITY simulation. Dotted lines are the 
volume integrals of the source terms on the right-hand side 
of Eqs [2] and ^] In steady-state, the solid and dotted lines 
should match. The small discrepancy near the outer edge of 
the simulation domain is likely due to numerical inaccuracy 
in flux calculations at the boundary. 



TABLE II: Analysis of TRINITY profile fits. S w and S Wl are 
the relative errors in total and stored energy, respectively. 
a is the RMS relative error associated with the subscripted 
quantity. 

Shot 5w <W 7 0~ n CTTi &T e 

JET 19649 0.09 0.12 0.11 0.12 0.13 
JET 42982 0.06 0.14 0.13 0.11 0.12 
AUG 19649 0.16 0.29 0.13 0.16 0.06 



profiles obtained in TRINITY, but not with the profiles 
taken from TRANSP. As an example of the sort of diag- 
nostic information available in our multiscale gyrokinetic 
calculations, we show the steady- state fluctuation am- 
plitudes for density, temperature, and electrostatic po- 
tential in Fig. [5] All fluctuation amplitudes are small 
compared to equilibrium quantities, giving us a check on 
our 5f <C / assumption. Interestingly, we see peaked 
electron temperature fluctuations at the magnetic axis, 
while the electrostatic potential increases with radius. 

Next, we consider JET shot #42982 1 EzEH wmc h 
achieved a record fusion energy yield of 22 MJ with 
21.5 MW of neutral beam heating. The plasma in this 
shot was a 50-50 deuterium-tritium mixture, operating 
in ELMy H-mode. Again, GS2 was used to calculate the 
turbulent fluxes, with a Miller local equilibrium model for 
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FIG. 5: Radial proflles for fluctuations of density, tempera- 
ture, and electrostatic potential calculated in GS2 for JET shot 
#19649. These fluctuations are obtained by time-averaging 
the instantaneous fluctuations computed in GS2 over the 
steady-state period at the end of the TRINITY simulation. 
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FIG. 6: Comparison of steady-state density and tempera- 
ture profiles constructed from JET shot #42982 by TRANSP 
(points and dotted lines) with those calculated in TRINITY 
(solid lines). 



the magnetic geometry. The perpendicular spatial grid 
used was the same as for the L-mode discharge, but the 
resolution in the remaining dimensions was increased to 
24 parallel grid points, 12 velocities, and 40 pitch angles. 
We again employed a hyperdiffusion operator to prevent 
a cascade to sub-grid scales. Each flux tube simulation 
was run 10 4 time steps, with average simulation times 
ranging from 250-1000 L^/vt^i- This time the electron 
density and ion and electron pressures were evaluated at 
10 radial locations, resulting in 40 flux tube simulations 
per transport time step. A total of 20 transport time 
steps were taken, and the full simulation ran just under 
ten hours on 8640 CRAY XT4 processors. 

As with the L-mode case, all profiles show relatively 
good quantitative agreement with their TRANSP counter- 
parts (Fig. [6|, with an RMS relative error averaged over 
all profiles of 12%. In this case, there appears to be a sys- 
tematic over-prediction of the ion and electron heat fluxes 
over the outer half of the minor radius, with the profiles 
pinned by the critical gradient (Fig. [7]). A possible ex- 
planation for this discrepancy is the fact that the radial 
flow shear, which has been found to have a significant 
effect in this discharge, 49 was not included i n our GS2 
simulations. This capability does exist in GS d 50 * 51 * and 
will be included in future TRINITY studies with evolving 
toroidal angular momentum. Again, fluctuation ampli- 
tudes are on the order of a percent of equilibrium am- 
plitudes, with electron temperature fluctuations peaking 
on-axis and electrostatic potential fluctuations increasing 
with radius (Fig. [8]). 

In order to study the effect of the edge temperature 
on the profiles, we repeated this simulation with electron 
and ion edge temperatures increased by 20%. The results 
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FIG. 7: Comparison of experimental (dotted lines) and simu- 
lated (solid lines) gradient scale lengths for JET shot #42982. 



are shown in Figs. 9ffTT We see that the 20% increase in 
edge temperature leads to an increase of approximately 
14% at the magnetic axis, as expected from the stiff pro- 
files indicated in Fig. [7[ The gradient scale lengths are 
similar to the base case across most of the minor radius, 
with the only significant discrepancies occurring near the 
edge where the stiffness of the profiles is less pronounced 
(Fig. 10). Fluctuation levels also do not differ signifi- 
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FIG. 8: Radial fluctuation profiles for density, temperature, 
and electrostatic potential calculated in GS2 for JET shot 
#42982. The fluctuation levels for this H-mode discharge 
are generally lower than for the L-mode discharge shown in 
Fig. [5] but they exhibit the same qualitative trends in their 
radial profiles. 



cantly from those obtained for the base case (Fig. 11). 

Finally, we consider ASDEX Upgrade shot #13151,^ 
which was an ELMy H-mode discharge with 5 MW of 
neutral beam heating. Here, we used GENE to calculate 
the turbulent fluxes, with a numerical magnetic equilib- 
rium generated by TRACER.^ We used 16 parallel grid 
points and a 64 x 48 grid in the perpendicular spatial 
plane, with perpendicular box widths of approximately 
64 pi at the outboard midplane. The velocity space was 
sampled with 32 parallel velocities and 8 magnetic mo- 
ments, resulting in a total of approximately 2.5 x 10 T 
mesh points for each GENE simulation. A linearized 
Landau-Boltzmann operator was used to model the ef- 
fect of collisions!^ The number of time steps taken in 
each flux tube simulation varied from about 3-4 xlO 4 de- 
pending on radial location, corresponding to simulation 
times of approximately 400-1000 L^Jv^^. Electron den- 
sity and ion and electron pressures were again evolved at 
8 radial locations, resulting in 32 flux tube calculations 
per transport time step. Within each transport time step, 
the total mesh points required was thus ~ 8 x 10 8 . Due 
to the increased spatial resolution in the perpendicular 
domain and the use of a physical collision operator, the 
GENE simulations took somewhat longer than the earlier 
GS2 simulations. In total, the simulation took just un- 
der 24 hours for 16 transport time steps on 16384 CRAY 
XT4 processors. 

The electron density and temperature agree relatively 
well across the minor radius, but the ion temperature is 
under-predicted near the magnetic axis (Fig. 12). This 
is reflected in the profile gradient scale lengths shown in 
Fig. 13 The RMS relative error averaged over profiles 
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FIG. 9: Comparison of steady-state density and temperature 
profiles for TRINITY simulations of JET shot #42982 with dif- 
ferent edge temperatures. The dotted lines with points corre- 
spond to a simulation using the edge temperatures reported 
by experiment, and the solid lines correspond to a simulation 
with the edge temperatures increased by 20%. We see that 
the increased edge temperatures lead to an increase in the 
temperatures near the magnetic axis of about 14%. 
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FIG. 10: Comparison of gradient scale lengths for two differ- 
ent TRINITY simulations of JET shot #42982. Dotted lines 
correspond to a simulation with edge temperatures taken from 
experiment, and solid lines correspond to a simulation with 
the edge temperatures increased by 20%. While the gradient 
scale length profiles are quite similar, the case with higher 
edge temperature leads to lower gradient scale lengths near 
the edge, where the profiles are likely less stiff. 
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FIG. 11: Radial fluctuation profiles for density, temperature, 
and electrostatic potential calculated in GS2 for JET shot 
#42982 with increased edge temperatures. They are quite 
similar to the case with edge temperatures taken from exper- 
iment. 
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FIG. 12: Comparison of steady- state density and tempera- 
ture profiles constructed from ASDEX Upgrade shot #13151 
by ASTRA (points and dotted lines) with those calculated in 
TRINITY (solid lines). 



is nonetheless only 12%. A lack of flow shear in the cal- 
culation of the turbulent fluxes is again a possible expla- 
nation for the discrepancy in ion temperature near- axis. 
Another possibility is the fact that no MHD model was 
used for the q profile computed by ASTRA, which resulted 
in q < 1 inside r/a ~ 0.4, with q(0) < 0.5. Consequently, 
there is significant uncertainty in modeling the system in 



CONCLUSIONS 



We have presented in this paper a complete theoret- 
ical and numerical model for the interaction of micro- 
and macro-physics in axisymmetric fusion devices. This 
model arises from a rigorous asymptotic expansion of 
the full Maxwell-Boltzmann system. The ordering as- 
sumptions used in the expansion were given in Sec. [TTJ 
along with the resulting closed system of equations for 
the evolution of the macroscopic thermodynamic profiles 
and magnetic geometry (Eqs. |2fl5l). These profiles de- 
pend on fluxes and heating (Eqs. ]5p| ) arising from classi- 
cal, neoclassical, and turbulent dynamics, requiring solu- 
tion of the neoclassical and gyrokinetic equations (Eqs. [9] 
and [lO ) . In order to evolve the magnetic geometry and 
close the system, the toroidal component of Ampere's 
Law (Eq. [l4|) and the Grad-Shafranov equation (Eq. 15) 
also have to be solved. 

In Sec. IIIIl we described the numerical scheme used 



in TRINITY to solve the multiscale gyrokinetic system 
of equations. The key idea in the approach is to use 
scale separations in space and time to embed a fine mesh 
for turbulent dynamics in a coarse grid for the transport 
solver. A local, nonlinear gyrokinetic code is used to cal- 
culate the turbulent fluxes, which are then passed to the 
transport solver. The macroscopic thermodynamic pro- 
files are then evolved in a manner consistent with the 
small p* limit of the orderings used to derive the non- 
linear gyrokinetic equation. Thus, TRINITY can be used 
to simulate time-dependent experimental phenomena. A 
Newton method is used to devise an implicit time step- 
ping algorithm, allowing for large time steps character- 
istic of the confinement time. With updated pressure 
profiles, one could then couple to solvers for the toroidal 
component of Faraday's Law and the Grad-Shafranov 
equation to evolve the magnetic geometry. 

Simulation results from TRINITY are provided in Sec. 
JV| with comparisons to JET and ASDEX Upgrade plas- 
mas. Relatively good agreement is found for all evolved 
profiles (density and electron/ion pressures), with an 
RMS relative error averaged over all profiles of 12%. 
Fluctuation levels for density, temperature, and electro- 
static potential are on the order of a percent across the 
minor radius, in general agreement with experimental ev- 
idence. 

While currently capable of faithfully simulating a range 
of interesting experimental conditions, there are several 
useful additions that could be made to the numerical 
model implemented in TRINITY. Coupling to a neoclassi- 
cal code which solves Eq. [9] would provide more accurate 
neoclassical fluxes, which can be important in certain ex- 
perimental regimes (for instance, when strong shear flow 
partially suppresses turbulent flux levels). Additionally, 
it would allow for the calculation of the parallel current 
necessary to evolve the safety factor and the poloidal flux. 
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FIG. 13: Comparison of experimental (dotted lines) and sim- 
ulated (solid lines) gradient scale lengths for ASDEX Upgrade 
shot #13151. 

Coupling to a Grad-Shafranov solver would then allow 
for a study of the effects associated with evolving mag- 
netic geometry. Additionally, coupling to a linear MHD 
code could be useful in monitoring macroscopic profiles 



to ensure that MHD stability boundaries are not crossed. 

In conclusion, we emphasize that the multiscale gy- 
rokinetic model presented here is not meant to be a com- 
prehensive model for all physics present in a tokamak 
discharge. Instead, it is meant to be used as a tool for 
studying the self-consistent interaction between micro- 
turbulence and macroscopic profiles. The ultimate goal 
of this approach is to obtain both a better qualitative and 
quantitative understanding of this interaction in order to 
enhance our ability to suppress turbulence and improve 
confinement in fusion experiments. 
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